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Abstract. The coalescence of massive black hole binaries is one of the main sources of low- 
frequency gravitational radiation that can be detected by LISA. When two galaxies containing 
massive black holes merge, a binary forms at the center of the new galaxy. We discuss the evolution 
of the binary after its separation decreases below one parsec. Whether or not stellar dynamical 
processes can drive the black holes to coalesce depends on the supply of stars that scatter against 
the binary. We discuss various mechanisms by which this supply can be replenished after the loss 
cone has been depleted. 

INTRODUCTION 

The prospect that low-frequency gravitational radiation will be detected by LISA has 
recently energized theoretical inquiries into the formation and the evolution of massive 
black hole binaries (MBHB). MBHB are sources of low-frequency gravitational radia- 
tion that will provide highest signal/noise for LISA, but the event rate for these sources 
is much less well known than for the other two principal astrophysical sources, com- 
pact binaries and extreme mass ratio inspiral (see contributions by G. Nelemans and 
S. Sigurdsson, this volume). Astronomical evidence for the existence of black holes 
with masses M\^> 10 6 M Q in galaxy spheroids with central stellar velocity dispersions 
o> 100 km s _1 is increasingly compelling; the evidence for black holes with masses 
100 — 10 M & in low-dispersion spheroids is still equivocal. When two galaxies merge, a 
MBHB forms at the center of the new galaxy QjL |2|] . There has been considerable inter- 
est in determining whether black hole coalescence occurs efficiently following galaxy 
mergers, since almost all predictions of MBH coalescence event rates equate the galaxy 
merger rate - derived from models of structure formation in which galaxies merge hi- 
erarchically iS 0, Bl 0. 0] - with the MBH coalescence rate. Detailed estimates of the 
MBH coalescence rate are fundamentally limited by the resolution of electromagnetic 
telescopes. It is therefore expected that MBH coalescences as observed by LISA will 
facilitate the first definite conclusions regarding galaxy merger rates and MBH demo- 
graphics at high redshifts [8]. 

The mutual gravitational capture of the black holes is facilitated by the dynamical 
drag imparted to the orbiting MBHs by the overlapping galaxies. The inner parts of 
elliptical galaxies and many spiral bulges are well described by power-law ("cuspy") 
stellar density profiles of the form p ~ r~ r with 7 ~ 2 BSD. High luminosity ellipticals 
often exhibit a shallowing of the density profile (7 < 2) near the very centers, where 
the MBHs are located. Numerical simulations of galaxy mergers have shown that the 
black holes remain embedded in their donor cusps throughout the merger II 1(11 . Since 



the dynamical drag is a function of the combined black hole and stellar cusp mass, 
which exceeds that of the black hole until the very final stages of the merger, the black 
hole inspiral and the MBHB formation take place on a dynamical time scale for black 
holes of comparable mass. In the unequal mass case, the infall time scale is lengthened 
proportionally to the black hole mass ratio M\jM<i (M\ > M 2 ), but for M 2 ^ 1O 6 M 
remains short compared to time scale for the subsequent evolution Il35|l . which is the 
subject of this contribution. 



GRAVITATIONAL SLINGSHOT AND BINARY DECAY 

The galaxy merger delivers the black holes to within the distance a hard = Gjl/Ao 2 , 
where }i = M\M2/{M\ + Mj) is the reduced mass. This distance is about 1 parsec 
for Mi = M2 = 2 x 1O 7 _M_0. At that point, the binary continues to decay by scattering 
stars super-elastically full . Stars in the merged galaxy with orbits approaching the 
binary closely enough to be perturbed by the rotating quadrupole component of the 
binary's potential belong to the "loss cone." The loss cone is defined in analogy with 
a similar structure characteri zing the distribution of the stellar-mass objects around 
solitary massive black holes [12, 13]. When a star inside the loss cone impinges on 
the MBHB, it exchanges kinetic energy with the binary and is shot out at an average 
velocity comparable to the binary's orbital velocity v* ~ a/ G{M\ +M2)/a ^> o. This is 
a form of the gravitational slingshot mechanism commonly used to accelerate spacecraft 
in the solar system. As a result, the binary's binding energy increases and its semi-major 
axis a decreases. 

In a crude approximation, the factor by which the binary separation decays can be 
related to the total mass in stars M scat that are scattered against the binary between times 
t\ and t2 via II 1411 : 

M scat (h,t 2 ) \ 



where ^ is the mass ejection coefficient which in fact depends weakly on a/ahard and 
the black hole masses. 

Approximation (HJ) is widely used in the semi-analytic modeling of the MBHB pop- 
ulation dynamics within hierarchical structure formation scenarios. Value of the mass 
ejection coefficient has been estimated from three-body scattering experiments 111411 : 
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J « 0.1 ln^— « 0.1-1. (2) 

This agrees with the value ^ ~ 0.5 measured in Af-body simulations of equal-mass 
MBH mergers II 1011 . If the galaxy is nearly spherical and the ejected stars return to 
interact with the binary more than once, the efficiency of mass ejection depends on 
the potential well depth A<J> separating the energies at which stars enter and exit the loss 
cone id: 

~ (A4>/2<j2)(a/a hard )- (3) 



To shrink the binary by one e-folding, a stellar mass of J 1 ' must be transported from 
an energy marginally bound to the black hole, to the galactic escape velocity. 



COALESCENCE 

If the semi-major axis decreases until it becomes less than lfl6ll : 

_ / 64 gMiMz(Mi +M 2 )F(e) \ 1/4 

«gr - | 5 c g % | , (4) 

the emission of gravity waves leads to the coalescence of the black holes on a time scale 
t gr . Here, F(e) is an eccentricity-dependent factor equal to unity for a circular binary. N- 
body simulations of MBHB formation in galaxy mergers suggest that the eccentricities 
remain moderate, although our understanding of the MBHB eccentricity evolution is 
still incomplete. The factor by which the binary must shrink from its conception to 
coalescence is: 



^hard 



where p = M2/M1 < 1 and a = dlogM^/dlogo ^4 — 5 is the logarithmic slope of 
the tight relation between the black hole mass and the central velocity dispersion of the 
galaxy OH- 

Therefore, to achieve coalescence in a Hubble time, the MBHB must shrink by 4 — 5 
e-foldings, which requires the scattering of 10 — 20 x ji worth of stars from the loss cone. 
After its initial emptying and the accompanying rapid shrinking of the binary by a factor 
of ~ 5 past ahard ( as observed in /V-body simulations with p ~ r~ 2 initial density profiles 
JkJ), the mass of the loss cone < /1 if the galaxy is approximately spherical. Once this 
mass is expended, the binary is still a factor > 10 wider than the separation at which it 
would coalesce gravitationally in a Hubble time; the binary decay may therefore stall 
I2III . The problem is exacerbated if the initial cusp profile is shallower 
than r~ 2 . The apparent inability of stellar dynamical processes to drive the binaries 
to coalescence poses a potential problem for the detection of these sources by LISA. 
Circumstantial evidence, however, suggests that long-lived MBHBs may be rare: no 
smoking-gun detections of MBHBs have been reported (with the possible exception 
of OJ 287; see W22I1 ). although a number of AGN have time-dependent features that 
have been provisionally attributed to MBHBs (Komossa, this volume). The theoretical 
difficulty of shrinking a MBHB by a factor of ~ 100 after its formation at a separation 
of ~ 1 pc is called the "final parsec problem." 



LOSS-CONE DIFFUSION 



The final parsec problem is most severe in nearly spherical galaxies where the mass 
inside the loss cone is the smallest. The loss cone boundary is defined by the minimum 



angular momentum 7i oss ~ yJG{M\ +M2)a that a star can have and still avoid being 
perturbed by the MBHB. Consider the state of the loss cone after all the stars inside 
have been scattered once. Since the depletion of the orbital population inside the loss 
cone leads to a stalling in the decay rate, continued decay hinges on the rate at which 
stars diffuse into the loss cone via collisional relaxation. The relaxation can be modeled 
by means of the orbit- averaged Fokker-Planck equation [23] describing the evolution 
of the phase space density f(E,J,t) subject to the boundary condition f(E,J\ OSSl t) = 0; 
i.e., stars are assumed to be removed from the system by the gravitational slingshot once 
they straddle the loss cone boundary. The loss cone flux is proportional to the gradient 
of the phase space distribution at its boundary, & oc df/ dJ\j loss . 

Although it is tempting to seek a steady-state solution /equi(^,^) 06 log(7/7i oss ) 11241 



25L 12111 . in reality the center of the galaxy is not likely to be collisionally relaxed ||2f 
and thus the distribution of stars near the loss cone is never in a steady state on time 
scales of order the relaxation time. Indeed, the final stages of a galaxy merger when the 
MBHB forms proceed in a time much shorter than the galaxy crossing time. Therefore 
the distribution function f(E,J) immediately following the formation of a hard binary 
can be far from that describing a steady-state flux of stars into the loss cone. Sudden 
draining of the loss cone during formation of the hard binary produces steep phase space 
gradients that are closer to the step function: 

f{E,J)n | fW> ■ ( 6 ) 

Since the collisional transport rate in phase space is proportional to the gradient of / with 
respect to 7, steep gradients imply an enhanced flux into the loss cone. The depletion 
of stars outside the loss cone affects the density profile of the galaxy and is identified 
with the cusp destruction. The broken power-law profiles of high-luminosity elliptical 
galaxies can be interpreted as fossil evidence for this process iEtLEsIi . 

The time evolution of the stellar distribution near the loss cone is an initial value 



problem equivalent to the diffusion of heat in cylindrical coordinates [12]. Ignoring the 



diffusion in E, the Fokker-Planck equation for diffusion in 7 reads (7 > J\ oss ): 

df(E,J,t) X(E) d 



dt J dJ 



\ df(E,J,t) 
dJ 



(7) 



where X(E) is related to the orbit- averaged diffusion coefficient. Since the boundary 
angular momentum 7i oss decreases with time, equation Q can be solved iteratively by 
discretizing the decrements in 7i oss and interpolating / between these via the Fourier- 
Bessel synthesis [15]. Solutions obtained this way can be compared to the collisionally 
relaxed, steady state (df/dt = 0) solution normalized to the isotropic distribution (Fig- 
ure^). In an example scaled to the galaxy M32 with a 3 x 1O 6 M black hole, the binary 
in the exact solution has decayed only ~ 30% more than that in the steady-state solution. 
The difference between the two, however, is much more substantial early on (Figure[IJ)), 
which is of crucial importance if episodic, violent dynamical perturbations such as the 
infall of satellite galaxies or giant molecular clouds rejuvenate the loss cone 112911 by 
restoring the steep phase-space gradients instrumental for the enhanced diffusion (Fig- 
ured]:). For example, if the episodic replenishment in a galaxy like M32 occurs every 




t (Myr) t (Myr) t 

FIGURE 1. (a) Evolution of the semi-major axis a in a merger of equal galaxies with ~ 10 6 M« black 
holes starting from the initial separation of a(0) = 0.1 pc (solid line). The evolution is always more rapid 
than that predicted assuming that the galaxy is collisionally relaxed (dashed line), (b) Enhancement of 
the decayed compared with that of the collisionally-relaxed, steady-state solution, expressed in terms of 
© = Aa/Afl equ i (triangles) and a power-law fit (solid line), (c) Schematic illustration of the evolution of 
the semi-major axis in the presence of episodic refilling of the loss cone. (From fl5ll .) 



10, 100, or 1,000 Myrs, the average MBHB decay rate will be « 10, 5, or 3 times 
higher than what the equilibrium theory would have implied 

Brownian motion of the MBHB in the neighborhood of the geometric center of the 
galaxy can to some extent mimic the effects of collisional relaxation and drive stars 
into the loss cone. The Brownian motion results from the equipartition of kinetic energy 
between the MBHB and the stars in the galaxy, (v^ rown ) ~ (m^/M^a 2 , where m* is the 
average stellar mass Il30l.l3 lL 13211 . As the binary wanders in space, it can sweep up stars 
that would have remained just outside the loss cone for a static binary. The time scale on 
which the loss cone refills in this fashion is: 

where K(E) is a function of the orbital energy such that K(2o 2 ) ~ 1. The amplitude of 
the Brownian motion is only modestly enhanced by "super-elastic scattering" by the 
binary ll30ll . In galaxies t^ mwn >l Gyr and thus the Brownian motion probably fails 
to significantly enhance the flux into the loss cone, in spite of earlier suggestions to 
the contrary IB3I1 . In Af-body simulations, however, the Brownian motion saturates the 
loss cone flux which is one of the many pitfalls that plague the numerical modeling of 
MBHBs. 



SPHERICAL AND ASPHERICAL GALAXIES 

Mode of interaction between the stars in a galaxy and a MBHB is also influenced by the 
geometry of the galactic potential, which could either be nearly spherical or substantially 
aspherical (axisymmetric, triaxial, or irregular). 

In spherical galaxies all stars that are candidates for slingshot ejection by the MBHB 
encounter the binary within one crossing time from the moment the binary forms. The 




FIGURE 2. (a) Snapshots of the binding energy distribution ^Y{E,t) of the stars residing inside the 
loss cone at t = 0. From right to left the data were taken at exponentially-increasing time intervals. As 
the binary separation decreases and the velocity of the slingshot ejection increases, stars inside the loss 
cone are heated through repeated scattering and shift to smaller E; a significant fraction of them remain 
inside the loss cone at all times, (b) Evolution of the semi-major axis exhibiting the logarithmic behavior 
motivated in the text (solid line). The slope of da/dlnt is close to that given by the theoretical prediction 
HOi : P(0) is the dynamical time of the potential well in which the galaxy is embedded (dashed line). 
(From CGI.) 



time-dependent loss cone solution derived in the previous section was based on the 
"sink" paradigm, in which a star is lost from the system as soon as it transgresses the loss 
cone boundary. This model is valid in the case of capture or tidal disruption of stars by a 
single black hole but is less relevant to MBHBs, since stars that interact with the binary 
simply receive kicks (AE,A7) that transport them to another orbit without necessarily 
ejecting them from the galaxy. If the binary orbit decays on a time scale longer than 
the orbital period of an interacting star, stars inside the loss cone can remain inside the 
loss cone after ejection, encountering the binary again at their next pericenter passage. 
In principle a star can interact many times with the binary before being ejected from the 
galaxy or falling outside the loss cone; each interaction takes additional energy from the 
binary and hastens its decay. 

We illustrate this "secondary slingshot" with a simple model in the spherical geom- 
etry. Consider a group of JV stars of mass m* and energy per unit mass E that interact 
with the binary and receive a mean energy increment of (AE). Averaged over a single 
orbital period P(E), the binary hardens at a rate II 1 511 : 

d (GM i M 2 \ (AE)jy 

j\-^r) =m *-pw (9) 

In subsequent interactions, the number of stars that remain inside the loss cone scales 
as 7j^ ss a, while the ejection energy is ~ Reject + G(M\ + Mij/la ~ cT l . Hence 
yy(AE) <x a l a~ l a . Assuming a singular isothermal sphere for the galaxy potential, 



we derive the result lfl5ll : 
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P(E )\ 

Hence the binary's binding energy increases as the logarithm of the time. Figure El 
illustrates the evolution of aT x in an /V-body simulation where star-star interactions 
have been replaced by a smooth potential to inhibit relaxation. The observed rate of 
decay da /dint ps 1.7 x 10 4 is close to the prediction of equation (flOb . where we 
have Ao 2 /G{M\ +M2) = 2 x 10 4 . Re-ejection at the rate of equation (fTOt would by 
itself induce changes in a by factors of a few over a Hubble time, in addition to the 
shrinkage due to collisional loss cone refilling. For P(E ) = 10 3 ' 5 ' 7 years, we obtain 
a(0)/a(?Hubbie) ~ 5, 4, and 3, respectively. 

In non- spherical galaxies where the total angular momentum J is not a conserved 
quantity, there exists a potentially larger population of orbits that can encounter the 
binary but only once per several orbital periods. The mechanisms of loss cone evolution 
and re-ejection discussed above in the spherical geometry would be modified somewhat 
in axisymmetric 12111 or triaxial galaxies [36]. The triaxial case is potentially the most 
interesting: stellar bars are commonly observed in galactic nuclei, and torques from 
barlike potentials are often invoked to drive gas inflows during the quasar epoch [47]. 
Since orbital angular momentum is not conserved in the triaxial geometry, a large 
fraction of the stars in a triaxial bar can potentially interact with the central MBHB. 
These "centrophilic" orbits are typically chaotic due to scattering off the central mass; 
in spite of their unfavorable time-averaged shapes, chaotic orbits can make up 50% or 
more of the total mass of a triaxial nucleus 13411 . 

Numerical integrations [36] reveal that the frequency of pericenter passages with 
r peri < a f° r a chaotic orbit of energy E is roughly linear in a, < a) ~ A(E)a, 

up to a maximum pericenter distance of r per i jmax (E) . The total rate at which stars pass 
within a distance a of the center is therefore: 

M^aJ A{E)^ ch&0 ,{E)dE (11) 

where ^ c \ mos (E)dE is the mass on chaotic orbits in the energy range E to E + dE. In a 
triaxial nucleus with density p ~ r~ 2 and central mass M^, the numerical integrations 
reveal: 

A(E)*5^.-^V, r^ S ^. (12) 
The feeding rate due to stars with energies E > 4>(r bh ) is then 

^. ° 3 a ~ _-l^chaos ( O V a 



M w ^chao s7 7— » 1O J M yr-'-j^ — — (13) 

G r bh 0.5 \200kms l J r hh 

with J^chaos me fraction of stars on chaotic orbits. Even a small J^chaos implies a 
substantial rate of supply to a MBHB when it first forms, with a ~ «hard ~ (M/2Mh>h)?"bh- 



Such high feeding rates would imply substantial changes in MBHB separations even 
if triaxial distortions to the potential were transient, due for instance to mergers or 
accretion events. 



Af-BODY SIMULATIONS 



Several attempts have been made to model the formation and the long-term evolution 
of MBHBs using large-scale A-body simulations MM MM. 

Because of the dis- 
creteness effects associated with approximating a galaxy with Af > 10 9 stars by a model 
consisting of only, at best, N < 10 6 particles, the applicability of direct A-body simula- 
tions appears to be limited to the early stage of the MBHB evolution. The rapidity with 
which a galactic merger proceeds guarantees that the discreteness effects do not affect 
the state of the galaxy just after the merger is complete. That point marks the transition 
to a much more gradual decay in the MBH separation when the effective relaxation time 
in the simulation can easily become shorter than the decay time. 

Simulations fail to reproduce the long-term evolution correctly because in the simula- 
tions the loss cone is almost completely full, while in real galaxies it is largely empty. In 
the general case there exists a critical stellar orbital energy E ct n separating the full region 
E < E cr - lt (large radii) from the empty region E > E CI n (small radii). Assuming a density 
profile p ~ r^ 2 and a potential of the form 2a 2 ln(r/ro) such that tq = IQ^GM^/la 2 
we find H: 



where J/^ = M^/m* is the number of stars that make the mass of the black hole. The 
transition from a full to an empty loss cone happens when E CI - lt becomes smaller than 
a few x2a 2 , implying that ~ 10 4 — 10 5 . Since a typical MBH contains 0.1% of 
its host galaxy's mass 13911 . and thus N ~ 10 3 =yfbh> an A-body simulation would have 
to contain 10 4 ~ 5 x 10 3 = 10 7 ~ 8 bodies to reproduce the correct, diffusive behavior of a 
real galaxy. This requirement is a severe one for direct-summation A-body codes, which 
rarely exceed particle numbers of ~ 10 6 even on parallel hardware 11321 . One route might 
be to couple the special purpose GRAPE hardware 1 , which is limited to Af <10 , to 
algorithms that can handle large particle numbers by swapping with a fast front end. 



We conclude by mentioning two outstanding problems related to the dynamics of MBHB 
binaries. Although the gravitational back-reaction tends to circularize binary black holes, 
some residual eccentricity can remain, especially if stellar dynamical processes prior 
to the emission of gravitational waves act to amplify the eccentricity to large values. 
Residual eccentricities result in the excitation of higher harmonics in the signal H40I1 . 



http ://grape . astron. s . u-toky o . ac j p/grape/ 




(14) 



OUTSTANDING PROBLEMS 



thereby complicating the detection of a gravitational wave event severely. In spite of 
the importance of eccentricity for detection, an accurate and general evaluation of the 
MBHB eccentricity evolution remains a challenge. 

Similarly, large-mass-ratio black hole binaries deserve extra attention. Although the 
best understood MBHBs are those of consisting of nearly equal-mass black holes, it is 
probable that most coalescences involve black holes of very unequal mass. For example, 
intermediate-mass black holes (IMBHs) with masses of lO 2 ~ 4 A/ may be able to form 
in young star clusters — such the Arches and the Quintuplet clusters in the Milky Way 
bulge 1141114211 — via the segregation of massive stars to the cluster center H43I1 . followed 
by the runaway growth in stellar mergers and the collapse of the agglomerate star into 
an IMBH II44L 14511 . An EVIBH makes its way to the galaxy center and forms a binary 
with the nuclear MBH 114611 . Therefore, large-mass-ratio MBHBs are expected to exist 
even in galaxies that had not experienced recent mergers. The orbital evolution of these 
systems remains another challenge to dynamical exploration. 



CONCLUSIONS 

We have focused on stellar dynamical mechanisms for extracting energy from a MBHB, 
but other schemes are possible and even likely. We note a close parallel between the 
"final parsec problem" and the problem of quasar fueling: in both cases, a quantity of 
mass of order ~ 1O 8 M must be supplied to the inner parsec in a time much shorter 
than the age of the universe. Nature clearly accomplishes this in the case of the quasars, 



probably through gas flows driven by torques from stellar bars [47]. The same inflow 
of gas could contribute to the decay of a MBHB, by leading to renewed star formation 
or rapid accretion of gas [48]. Similarly, the presence of a third black hole in a galactic 
nucleus could accelerate the decay by increasing the MBHB's eccentricity through the 
Kozai mechanism 114911. or by extracting the binary's energy via the triple black hole 
slingshot interaction IL50I1 . 
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